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In this article we investigate the hitting time of some given bound- 
aries for Bessel processes. The main motivation is coming from math- 
ematical finance when dealing with volatility models but the results 
can also be used in optimal control problems. The aim is here to con- 
^ ■ struct a new and efficient algorithm in order to approach this hitting 

time. As an application we will consider the hitting time of a given 



^ ( , level for the Cox-IngersoU-Ross process. The main tools we use are 



on one side an adaptation of the method of images to this particu- 
lar situation and on the other side the connection existing between 
Cox-IngersoU-Ross processes and Bessel processes. 



Oh ' 1. Introduction. The aim of this paper is to study the hitting time of some curved boundaries 

. for the Bessel process. Our main motivations come from mathematical finance, optimal control and 

, neurosciences. In finance Cox-Ingersoll-Ross processes are widely used to model interest rates. As 

an application, in this article we will consider the simulation of the first hitting time of a given 
level for the CIR by using its relation with the Bessel process. In neurosciences the firing time of a 
neuron is usually modelled as the hitting time of a stochastic process associated with the membrane 
\ potential behaviour (for introduction of noise in neuron systems see Part I Chapter 5 in [7]). The 

literature proposes different continuous stochastic models like for instance the family of integrate- 
. and- fire models (see Chapter 10 in [6]). Most of them are related to the Ornstein-Uhlenbeck process 

I which appears in a natural way as extension of Stein's model, a classical discrete model. In Feller's 

model, generalized Bessel processes appear as a more realistic alternative to the Ornstein-Uhlenbeck 
process, see for instance [10] for a comparison of these models. Therefore the interspike interval, 
which is interpreted as the first passage time of the membrane potential through a given threshold 
is closely related to the first hitting time of a curved boundary for some Bessel processes. 

Our main results and the main algorithm are obtained for the case of Bessel processes. We use 
in our numerical method the particular formula that we obtain for the hitting time of some curved 
I boundaries for the Bessel process and the connection that exists between a Bessel process and the 

Euclidean norm of a Brownian motion when calculating the hitting position. As an application 
we consider the hitting time of a given level for the Cox- Ingersoll- Ross process. In order to obtain 
this, we use first of all the connection that exists between CIR processes and Bessel processes and 
secondly the method of images for this particular situation. 

The study of Bessel processes and their hitting times occupies a huge mathematical literature. 
Let us only mention few of them: A. Going-Jaeschke and M. Yor [8] consider a particular case 
of CIR processes which are connected with radial Ornstein-Uhlenbeck processes and their hitting 
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times; L. Alili and P. Patie [1] investigate as a special situation the Bessel processes via some 
boundary crossing identities for diffusions having the time inversion property; recently P. Salminen 
and M. Yor consider the hitting time of affine boundaries for the 3-dimensional Bessel process [16]. 

In a recent paper Y. Hamana and H. Matsumoto [9] give explicit expressions for the distribution 
functions and the densities of the first hitting time of a given level for the Bessel process. Their 
results cover all the cases. Let us also mention a recent work of T. Byczkowski, J. Malecki and M. 
Ryznar [2] . By using an integral formula for the density of the first hitting time of a given level of 
the Bessel process, they are able to obtain uniform estimates of the hitting time density function. 

In all these papers the formulas are explicit and are hard to use for a numerical purpose as 
they exhibit Bessel functions. The main idea of this paper is to get rid of this difficulty by using 
two important tools: first of all the method of images that allows to obtain, for some particular 
boundaries, an explicit form for the density of the hitting time and secondly the connection between 
5-dimensional Bessel processes and the Euclidean norm of a (5-dimensional Brownian motion in order 
to get the simulated exit position. By coupling these ingredients we are able to construct a numerical 
algorithm easy to implement and very efficient which approaches the desired hitting time. 

We will use here a modified version of the random walk on spheres method which was first intro- 
duced by Muller [13] in 1956. This procedure allows to solve a Dirichlet boundary value problem. 
The idea is to simulate iteratively, for the Brownian motion, the exit position from the largest 
sphere included in the domain and centered in the starting point. This exit position becomes the 
new starting point and the procedure goes on until the exit point is close enough to the boundary. 
Let us notice that the simulation of the exit time from a sphere is numerically costly. 

The method of images was introduced in 1969 by H.E. Daniels [4] as a tool to construct nonlinear 
boundaries for which explicit calculations for the exit distribution for the Brownian motion are 
possible. The method was developed also in H. R. Lerche [11]. We adapt this method for the Bessel 
process by using the explicit form of its density. For some particular curved boundaries we can 
explicitly evaluate the density of the Bessel hitting time. 

The paper is organized as follows. First we present some new results on hitting times for Bessel 
processes. Secondly we construct the new algorithm for approaching the hitting time, the so called 
Walk on Moving Spheres algorithm. Finally we present some numerical results and as a particular 
application the evaluation of the hitting time of a given level for the Cox-Ingersoll-Ross process. 

2. Hitting time for Bessel processes. 

Bessel processes play an important role both in the study of stochastic processes like Brownian 
motion and in various theoretical and practical applications as for example in finance. 

Let us consider the 5-dimensional Bessel process starting from y, solution of the following stochas- 
tic differential equation: 



the index of this process. We call 6 the dimension of the process. This terminology is coming 
from the fact that, in the case of positive integer 5 G N, a (5-dimensional Bessel process can be 



(2.1) 




where {Bt)t>o is a one dimensional Brownian motion. We denote 



(2.2) 
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represented as the Euclidean norm of a 5-dimensional Brownian motion. This will be a key point 
in our numerical method. 

The density of this process starting from y is given by: 

(2.3) py{t,x) = '^{^ exp (^-^^^^ (^) , fort>0, y>0,x>0, 
where Iu{z) is the Bessel function whose expression writes: 

oo 

(2.4) ^-(^) = E©' 



n=0 

When starting from y = 0, the density of Z^'^ is: 



n!r(i/ + n + 1) 



11 1 f 

(2.5) Po{t,x) = --— -x'^-^exp , for t > 0, X > 0. 

2.1. The method of images for Bessel processes. In this section, we investigate the first hitting 
time of a curved boundary for the Bessel process starting from the origin. Let ^(t) denote the 
boundary, and introduce the following stopping time: 

(2.6) = inf{t > 0; ° > ij{t)}. 

For some suitable choice of the boundary, the distribution of can be explicitly computed. The 
idea is given by the following remark on the method of images (see for instance [4] for the origin of 
this method and [11] for a complete presentation): 

Fundamental idea. Suppose that F is a positive cj-finite measure satisfying some integrability 
conditions (to be specified later on) and define 

(2.7) u{t,x) = po{t,x) - - I py{t,x)F{dy), 

for some real constant a > 0. Let 

tp{t) = inf{x G M; u{t, x) < 0}, for ah t > 0. 
Then u{t, x) is solution of the partial differential equation: 



(2.8) 



■u(t,V'(i)) = 0, foralH>0, 
[ u{0,.)=6o{.) on (-00,^(0+)]. 



From this remark we deduce an interesting expression for the hitting time. We can prove that: 

= inf{t > 0; u{t, zf'°) = 0}. 

This means simply that in order to obtain informations on the hitting time it sufficies to look for 
u{t, °) = 0. 

Let us express this in a general result 



Theorem 2.1. Let F{dy) be a positive a-finite measure such that Po(t, ^/ey)F{dy) < oo for 
all e > 0. Let a > and define the function: 



(2.9) 



u{t,x) = po{t,x) / py{t,x)F{dy). 



Consider il){t) such that u{t,tp{t)) = 0. Then the probability density function of is given by 

5-1 



(2.10) 



\{t^ G dt) 



Idu , 



x=4){t) 



Idu , 



x=0 



2x 



-u{t, x) 



x=0 



dt. 



Proof. We will only point out the main ideas for the proof in this case as it follows mainly the 
ideas introduced in [11]. A complete description of the method and this result for the Brownian 
motion case can be found in [11]. 

Let us consider 

(2.11) uit, x) = po{t, x)-- [ py{t, x)Fidy), 

where F{dy) is a measure on We consider tp{t) the solution of u{t,tp{t)) = 0. Let us define as 
before = inf{t > 0; zf'^ > ip{t)}. Then u{t,x)dx = P(r^ > t, G dx) and 



(2.12) 



^{t^ > t) = / u{t,x)dx. 
Jo 



In order to get the distribution of we have to evaluate the derivative of Po(t^ > t). By using the 
equality (2.12) we obtain: 



"o(t^ G dt) 



(2.13) 



-i;'{t)u{t,m) 

1 /-^w^^n 



m du 



{t, x)dx dt 



dx^ 



{t, x)dx 



dt 

6 - 1 /■'^(*) d_ ( 1 
dx 



u{t, x) I da; dt, 



as u{t,x) is solution of the partial differential equation (2.8). We obtain thus: 
(2.14) 



Mr^ e dt) 



19m 
2dx 



{t,x] 



1 du 
2dx 



{t,x) 



6-1 
2x 



u(t, x) 



dt. 



x=0, 



As 2^(l-j u{t, '4^{t)) = and this ends the proof of the theorem. 



□ 



The idea behind the method of images is that for some particular forms of F{dy) we can derive 
explicit formulas of the hitting time distribution. More precisely: 



Proposition 2.2. Let, for 5 = 2u + 2> Q and a > 
2.15) 



2Hog: 



V-aW - Y""'^^^r(z. + l)t-+i2- 
Then the probability density of is given by 



for t < 



_r(i/ + 1)2^ 



1 



(2.16) 



1 



^^^^ = ^t{ 2* log 



T{i^ + l)t''+^2'- 



dt. 



Proof. By using the expression in (2.3) we remark first that 

(2.17) y2.+i^,(t,x) = x2'^+V.(t,2/). 

Let us consider, as in Theorem 2.1 

1 



(2.18) 



u{t,x) =pQ{t,x) 



Py{t,x)F{dy), 



with F{dy) = y ~^ l^y~^Qydy. In this situation the function u defined in (2.18) writes 



1 



(2.19) 



u{t, x) = po{t, x) X 

11 1 



2u+l 



■ exp 



_ 1 

2t / ~ a ' ^ 



2u+l 



2" t'^+ir(i/ + i) 

For simpHcity we wih write instead of tpa- Following the result in the Theorem 2.1, we are looking 
for ip{t) such that u{t,ip{t)) = 0. This yields : 



(2.20) 



X = tpit) = . 2t\og 



under the obvious condition t^'^ < 
We can now notice that: 



Po{t,ilj{t)) = -(V(i)) 
a 



2u+l 



and we can prove easily that: 



pit,x) = i6-l)^-^-Mt,x). 
ox X t 



We obtain, after replacing in (2.10) and after applying the Theorem 2.1, for this particular case: 

Foir^ edt) =j ^P{t)po{t, ^P{t))dt 



T{iy + l)t^+i2'^ 



dt. 



which gives the desired result. 



□ 



The second boundary which allows to express explicit results is obtained by using the Markov 
property for the Bessel process. 

Proposition 2.3. Let, for 5 = 2z/ + 2>0, s>0 and a > fixed, 
(2.21) 



2t{t + s) 



(z^ + l)log (1 + - ) + loga 



for all t > if a > 1 and for t < ^ if < a < 1. Then the probability density function of 

the hitting time is given by: 
(2.22) 



1 I ft + sY 



log a 



t + s 



u+1 



exp 



t + s 



log a 



t + s 



dt. 



Proof. We will only sketch the proof as it follows the same ideas as the one of the Theorem 
2.1. Let us consider the measure F{dy) = pQ{s,y)dy for s > fixed. Then, when evaluating the 
corresponding u{t, x) we have: 



u{t,x) =pQ{t,x) 



Po{s,y)py{t,x)dy 



Poit,x) po{t + s,x) 

a 



1 



1 



^21/+! 



expl-- 



1 



1 



a{t + s)" 



+1 



exp 



2{t + s) 



by using the Markov property. We obtain the form of '(/'(i) by the condition u(t, = which 
gives: 



(2.23) 



2t{t + s) 



(i/ + 1) log ( 1 + - ) + log a 



for t > if a > 1 



for t < 



1 if a < 1. 

(1)PTT_1 



In order to obtain the distribution of one has only to evaluate: 



(2.24) 



dx^ ' ' ^ ' x t{t + s) 



xpo{t,x) 



and ^i&^il for x = and x = 'il)[t) and replace the values in the general form (2.14). The expression 
(2.22) follows. □ 

Remark 2.4. We can notice that the function ipait) defined by (2.21) satisfies for large times 

i^a{t) ^ y/t fora=l, 
il^ait) ^ t for all a > I. 

In particular, we can approach large times by considering this kind of boundary. 

A new boundary can be obtained by using the Laplace transform of the square of the 6 - 
dimensional Bessel process starting from x. More precisely: 

Proposition 2.5. Let, for 6 = 2u + 2 > 0, A>0 and a > fixed, 

Xt 



(2.25) 

for 

(2.26) 



l + 2Xt 



t\ 



A 



2 , a(l + 2Xtr+^ 
+ - log 



i + 2Aty t ^ 2'^t^+ir(i/ + 1) 



t < 



1 /2^r(z^ + 1)\~ 



t > 



^/A> 



1 /2''r(z^ + 1)\^ 



Then the probability density function of the hitting time is given by: 
(2.27) Po(t^ G dt) = 1 



A 



2^ a{l + 2XtY+^ , ,,,,, 



l + 2XtJ t 2^t^+ir(i/ + 1) 
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Proof. We present only the main ideas as the resuh follows as above from the general method 
in Theorem 2.1 applied to the measure F{dy) = y^i^+ie""^^ ^{y>o}'^y- For this measure u{t,x) takes 
the form: 

- I pyit,x)F{dy) 
1 



u{t,x) = po{t,x) 
= Po{t,x) 



Py{t,x)y 



dy 



Po{t,x) 



-E e 



By using the expression of the Laplace transform for Z^'^ we obtain 



(2.28) 



u{t,x) = poit,x) 



1 



exp 



Ax 



l + 2Xt 



a (1 + 2Xty/^ 

We consider first the equality u{t,^lJ{t)) = in (2.28) and this gives the form of ^l^{t) in (2.25). 
Afterwards, we can evaluate once again in this particular situation 

Xt 



^{t,x) = i6-l)^ 
ox X 



po{t,x). 



t l + 2Xt^ 

For this particular case, there is only one no vanishing term in the expression (2.10) of Po(t^ G dt), 
that is the term — (^j — r^^At) Po{t^ ^) of |^(*> ^) ^ = this is exactly given by the right 

side of the formula (2.27). □ 

Corollary 2.6. The previous results writes for 6 = 2 



(1) For a > 0, t < a and '4'{t) = ^/2tlog — , the density of the hitting time is 

Po(T^edt) = ^log^dt. 

2a t 



(2) For s > 0, a > 0, t < ^ and tpit) 
function of is given by: 



'2t(t + s) 



log a 



t + s 



, the probability density 



Po(r^ G dt) 



(3) For a > and ^(t) 



(2.29) 



At 



1 + 2At 



t + s 



+ t^ 



log a 



t + s 



exp 



t + s 
t 



log a 



t + s 
1~ 



dt. 



X 



l + 2Xt 
t > 0, 



2, a(l + 2At) ^ 
+ ylog ^ ' for 



1 



ifX>—, 
■' - 2a 



t< ^ i/A<— , 

- 1 - 2Ao 2a 



the probability density function of is: 
Po(r^ G dt) = 1 



A 



2, a(l + 2At) , ,,,,, 
l + 2At/ +7log^^Po(t,^(t))dt. 
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2.2. Approximation of the first hitting time for Bessel processes starting from the origin. In 
this section we will construct a stepwise procedure, the so-called random Walk on Moving Spheres 
(WoMS) algorithm, which allows to approach the first time the standard Bessel process hits a 
given level I > 0. Of course, this stopping time ti = inf{t > 0; Z^^ = 1} can be characterized by its 
well-known Laplace transform computed by solving an eigenvalue problem. Indeed if {Z^'^ , t > 0) 
is the Bessel process starting from x, of index ly = i — 1 then, for > and x <l, we get 



e 



Here 1^ denotes the modified Bessel function. This Laplace transform can be used to describe the 
tail distribution: Ciesielski and Taylor [3] proved that, for (5 = 2i/ + 2 G N*, 

ff^o T« >t)= > - — e ^ . 

2^" ^r(z^ + ^) ^ X+i(>,fe) 

where J7i/ is the Bessel function of the first kind and {ju,k)u,k is the associated sequence of its positive 
zeros. 

We are looking for a numerical approach for the hitting time and these formulas are not easy 
to handle and approach, in particular we can not compute the inverse cumulative function ! The 
aim of this section is to construct an easy and efficient algorithm without need of inverting the 
Laplace transform and without using discretization schemes since the hitting times are unbounded. 
In the next step we will extend this procedure to the hitting time of time-dependent boundaries 
like straight lines, useful in the description of the hitting time of a given level for the CIR process 
(the Laplace transform is then unknown). 



2.2.1. Hitting time of a given level for the Bessel process with positive integer dimension 5. Let 

(k) 

US consider 6 independent one-dimensional Brownian motions [B^ ,t'>0),l<k<6. Then the 
Bessel process of index 6 starting from 0, satisfies the following property: 




{Zp'', t > 0) has the same distribution as 
Let 

(2.30) Ti = inf{t > 0; ° > /}. 

In particular, we can express Ti by using the first time when the 5-dimensional Brownian motion B = 
{B^^\ . . . , B^^'^) exits from the Euclidean ball D centered in the origin with radius /. Approximating 
the exit time and the exit position for the 2-dimensional Brownian motion of a convex domain was 
already studied by Milstein [12]. He used the so-called random walk on spheres algorithm which 
allows to approach the exit location and the exit time through an efficient algorithm. The exit 
position is really simple to obtain (as it is uniformly distributed on the circle) while the exit time 
is much more difficult to approach. That is why we will construct an adaptation of this initial 
procedure in order to obtain nice and efficient results concerning the Bessel process exit time. 

Let us introduce now our Walk on Moving Spheres (WoMS) algorithm. We first define a continuous 
function p -.M? ^ which represents the distance to the boundary of D: 

(2.31) p{x) = inf{||x - y\\; y e D"} = I - \\x\\. 



For any small enough parameter e > 0, we will denote by the sphere centered at the origin with 
radius I — e 

(2.32) D' = {x£ D; \\x\\ < I - e} = {x £ D; p{x) > e}. 



Algorithm (Al) for 5 = 2. Let us fix a parameter < 7 < 1. 
Initialization: Set X(0) = (Xi(0), X2(0)) = 0, 6^0 = 0, Go = 0, ^0 = 

First step: Let (Ui, Vi, Wi) be a vector of three independent random variables uniformly distributed 
on [0, 1]. We set 

h = AoUiVi, 91 = 90 + ^1, 

' cos(27rT^i) " 



x{iy = (Xi(i),X2(i))T = x{oy + VAo(ei) 



sm{2TTWi) 



where 



(2.33) ^pa{t) = Y'2tlog ^, t < o, a > 0. 

At the end of this step we set Ai = j'^ p(X{l))'^e/2. 

The n-th step: While X[n — 1) G , simulate {Un, Vn, Wn) a vector of three independent random 
variables uniformly distributed on [0, 1] and define 



(2.34) 



On — An-lUnVn, 9.„ — 9„_i + 9n, 

X{7iy = {Xi{n),X2{n)y = X{n - 1)t + V^„_,(0„) 



cos(27rW„ 
sin(27rW„ 



At the end of this step we set A^ = 7^p(X(n))^e/2. 

When X{n — 1) ^ the algorithm is stopped: we set On = 0, 9„ = 9„_i and X[n) = X{n — 1). 
Outcome: The hitting time 9„ and the exit position X{n). 



Remark 2.7. The WoMS algorithm describes a D-valued Markov chain (X(n), n > 0). Each 
step corresponds to an exit problem for the 2-dimensional Brownian motion. If X{n) = x then we 
focus our attention to the exit problem of the ball centered in x and of radius (ipait), t >0): the exit 
location corresponds to X{n + 1) and the exit time to 6n+i- Of course the choice of the parameter a 
plays a crucial role since the moving sphere has to belong to the domain D as time elapses... When 
the Markov chain X is close to the boundary dD, we stop the algorithm and obtain therefore a good 
approximation of the exit problem ofD. 

Comparison with the classical (WoS) algorithm: at each step, the n-th step of the classical walk 
on spheres (W oS ) is based on the exit location and exit time, which are mutually independent, for 
the Brownian paths exiting from a ball centered in X{n— 1) and with radius ^p{X{n— 1)). The exit 
location is uniformly distributed on the sphere while the exit time is characterized by its Laplace 
transform. Therefore, if one knows X(n — 1) then the diameter of the sphere is deterministic. For 
the WoMS the center of the ball should also be X{n — 1) but the radius is random, smaller than 
^p[X{n — 1)). The exit location will also be uniformly distributed on the sphere but the exit time 
will be much easier to simulate: in particular, you don't need to evaluate the Bessel functions. 

9 



The stochastic process {X{n), n > 0) is an homogeneous Markov chain stopped at the first time 
it exits from D^. In the following, we shall denote this stopping time which represents in fact 
the number of steps in the algorithm: 

iV^ = inf{n > 0; X{n) ^ D^}. 

We just notice that X{N^) £ by its definition. 

The algorithm (Al) is presented in the 2-dimensional case. Of course we can construct a generaliza- 
tion of this procedure for the (5-dimensional Bessel process when 5 € N*. For notational simplicity, 
we use a slightly different method: instead of dealing with a Markov chain (X(n), n G N) living in 
R'' we shall consider its squared norm, which is also (surprisingly) a Markov chain. At each step, 
we shall construct a couple of random variables xi^)) associated to an exit problem, the first 
coordinate corresponds to an exit time and the second one to the norm of the exit location. 
We introduce some notations: represents the unit ball in M'', and vri : — t- M the projection on 
the first coordinate. 



Algorithm (A2). Let us fix a parameter < 7 < 1. 
Initialization: Set x(0) = 0, = 0, Eq = 0, ^0 = (7^^^e/(i/ + 

The n-th step: While \J x{n — 1) < I — e, we choose Un a uniform distributed random vector on 
[0, Gn a standard gaussian random variable and Vn an uniformly distributed random vector 

on . Consider Un, Gn and Vn independent. We set 



(2.35) 

where 
(2.36) 



= {r^TiW ^"(1) • • • ^"(L^J + 2)) exp { - }, 

^ Xin) = x{n - 1) + 27ri(K)Vx(^-l)V'A_i(e«) + V^l.^lCn), 



"n— 1 ~l~ S,n, 



2tlog 



r(i/ + i)t'^+i2'' 



t < tmax(a) : 



At the end of this step we set 

An=(7\l-Vx{n)fe/{,^ + l) 



r(i/ + 1)2^ 



1 



, a > 0. 



When sj x("') > ^ — £ ih.^ algorithm is stopped: we then set ^„ = 0, H„ = and xin) = xin — !)• 

Outcome: The hitting time H„ and the value of the Markov chain x(^)- 



It is obvious that for the particular dimension 6 = 2, that is = 0, the stopping times obtained by 
the algorithm (Al) and (A2) have the same distribution. Moreover, for each n, x(n) has the same 
distribution as ||X(n)|p. In other words, the number of steps of (Al) and (A2) are identical in law, 
the number of steps will be denoted in both cases A'^^. 

Theorem 2.8. Set 6 £ W. The number of steps of the algorithm WoMS (A2) is almost 
surely finite. Moreover, there exist a constant > and £q{5) > 0, such that 

HN^] < Gs\ logel, for all e < eo{6). 
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Theorem 2.9. Set 6 & W . As e goes to zero, Hate converges in distribution towards ti the 
hitting time of the 6-dimensional Bessel process (with cumulative distribution function F), which 
is almost surely finite. Moreover, for any a > small enough, 

(2.37) f 1 - -j^)F%t -a)< F{t) < F%t), for all t > 0, 



\/2mr 

where F^{t) := F{En^ < t). 

These results and the key ideas of the proofs are adapted from the classical random walk on 
spheres (WoS), see [12]. 
Proof of Theorem 2.8. 

Step 1. Let us estimate the number of steps. Since (x(^)) ?^ > 0) is a homogeneous Markov chain, 
we introduce the operator Pxf defined, for any non-negative function / : — t- M+, by 

P.f--= I /(y)P(x,dy), 

where P(x, dy) is the transition probability of the Markov chain. By definition, xi'^ + 1) depends 
only on xip)-, and Let us note that, by construction, Vn and ^„ are independent. Moreover 

using the result developed in the Appendix, the density of ^n{ ^ a^'^^^^ ) ^® given by: 

(2.38) ii{r) = ^ogr) l[o,i](r). 

If we denote a'^ the uniform surface measure on the unit sphere in M*^, we get 

(2.39) p^f = j f(x + 2TTi{u),/^K{x,r) + K'^{x,r)^fi{r)draf{du), 
with K{x,r) defined by 

(2.40) K{x,r) = iPA^^ ^ 



1 

r 



and A depending on x in the following way: 

'72(/-^)2e\<^+ir(i/ + l) 



A 



We can observe the following scaling property: V'A(^'^+^i) = A'^''+'^ 'il)i{t) . Therefore the definition 
of ipi leads to 



(2.41) K 



(x,r) = 7(/ - V^)W-— pylog-^ = 7(Z - ^/^)^/er{-\ogr). 



Step 2. Using classical potential theory for discrete time Markov chains (see, for instance, Theorem 
4.2.3 in [14]), we know that 

cP{x) =^Jy, 9{x{n))\ 

\ n=0 / 
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satisfies, for any non-negative function g, 

j (l)ix) = P^(l) + gix), 0<x<(Z-e)2, 
^ ^ ^ I <P{{l-er) = 0. 

In particular, for g = 1, we obtain that 4>{x) = Ex[N^]. In order to get an upper-bound for the 
averaged number of steps, it suffices to apply a comparison result. We choose the function 

(2.43) U%x) = {log((/ - ^)/e) - log(l - 7)}/{Cs7^) 

which satisfies U'^{x) > PxU^ + 1, for all < x < (/ — e)^ (see Lemma 2.10 for the definition of the 
constant and for the inequality) and U^{x) > for all < x < (/ — e)^. A classical comparison result 
related to the potential theory (see, for instance, Theorem 4.2.3 in [14]) implies that E2;[A^^] < U^{x) 
for all X E [0, {I — e)^] and consequently leads to the announced statement. □ 

Lemma 2.10. Let us define, for small e > 0, C/^(x) = {log((/ - \/x)/e) - log(l - j)}/{Csj'^) 
for X G [0,/^[, where 7 is related to the definition of the WoMS. There exists a constant Cg > 0, 
such that, for any x g]0, (/ — e)^[, the following inequality yields 

PxU' - U'{x) < -1. 

We recall that P^U^ is defined by (2.39) and (2.41). 

Proof. We will split the proof in several steps. 

Step 1. First of all, we observe that > — log(l — 7)/(C57^) in the domain [0, (Z — e)^]. Let us 
consider now x(0) = x S [0, (/ — e)^] and y in the support of the law of x(l) and let us prove that 
U^{y) > 0. By the definition of x(l) we obtain 

X(l)< sup {x + 2y^i,A{t)+i^\{t)), 

max 

where ^= (^"^(l - ^fe/ {u + l)Y and both and t 

max are defined by (2.36). The right 
hand side of the preceding inequality is increasing with respect to y so that: 

X(1)<(V^+ sup ^A{t)'^ 

i6[0,imax{A)] 



Furthermore, for a > the maximum of the function is reached for t^ 
and is equal to 



1 

1 / a \ "+1 



e \ T{u+l)2^ 



1 >. 1/2 

2{u + 1) ( a \ -+i I 



(2-44) sup Va(t) , \v( 

Finally using the definition of A and the inequality x < {I — e)^, we find the following lowerbound: 



I - VxO) > (^ - ^/^)(l - 7) > ^(1 - 7)- 

We can therefore conclude that, for any y in the support of the law of x{V) (even for y > {I — e)^), 
U^{y) > which ensures that is well defined and non-negative in the domain of the operator 
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Step 2. Furthermore the Taylor expansion yields 

(2.45) U{y)<U (X) + _ - _ + ^^^^^^^ _ x,yG[0,l[. 

If x(0) = 2; and y is in the support of the random variable x{^)j then 



^/y — \fx = Y X + 2'K\[u)\fxK(x^ r) + i^^i^^;^ j.'j _ 
> 7ri(ti)A'(x, r). 

By the expansion (2.45) and the definition of the operator Px given by (2.39), the following upper- 
bound for the operator Px holds 

P^U' = j (x + 2Tri{u)y/^K{x,r) + K^{x,r)^ fi{r)dr af{du), 



< W{x) - f f 
Jo Js 



Tri{u)K{x,r) . 



jsi 2Csi\i - v^y 

1 f ■K'f{u)K^{x,r'' 



//(r)dr ai (du) 

lo Js' ^CsT{l - Vx)-^ 
where 

(2.46) si := {ueS^ : Tri{u) > 0}. 

Due to symmetry properties, the first and the third integral terms vanish. Then (2.41) leads to 

P.U' <U'{x)-^ [ 7Tl{u)af{du) 

with 

(z^+ir+2e /■V-+i(_logr)-+2dr. 

The description of the probability density function in the Appendix leads to the following explicit 
value 



v + 1\ e 



^iy + 2j T{u + 2) 
In order to conclude the proof, it suffices to choose the particular constant 

Cs = I [ ttUu) afidu). 
Jsi 



□ 
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Proof of Theorem 2. 9 

The proof is sphted in two parts. First the steps of the algorithm and the hitting time of the 
Bessel process of index i/ shall be related to stopping times of a (^-dimensional Brownian motion 
(z^ = I — 1). Secondly we point out that the corresponding stopping times are close together by 
evaluating deviations of the Brownian paths. 

Step 1. Let B = {B^^\ B^'^\ . . . , B^^^) be a (5-dimensional Brownian motion. Then the norm of B 
has the same distribution as a Bessel process of index (see, for instance [15]). Hence the first 
hitting time r; is identical in law to the stopping time 

Ti = mf{t > 0; Bt ^ D}, 

where D is the Euclidean ball centered at the origin and of radius I. We introduce then a procedure 
in order to come close to T;. For the first step we shall focus our attention to the first exit time 
of a moving sphere centered at the origin and of radius ipa{t) defined by (2.36), we denote this 
stopping time. Of course this moving sphere should always stay in D, so we choose a such that the 
maximum of ipa stays smaller than /. By (2.44), we get 

, m , r(z. + i) / eP 

sup -0^(0 < I a < — — 
t<a z + 

For a = Aq = ^^^2^^ i^ V+i ) with a parameter 7 < 1, the condition is satisfied: sup^^ Va(0 = 

^2u+2i ^ describe the law of (^i,B^-^). The norm of the Brownian motion is identical in 

law with the Bessel process; therefore Proposition 2.2 implies that the density function of .^i is 
given by (2.16) with a replaced by ^o- Using the law described in Proposition A.l, we can prove 
that ^1 has the same distribution as 

An \ ^ „z 
e , 



r(z^ + l)2 

where Z is Gamma distributed with parameters a = iy + 2 and p = ^^q-j- . By construction we deduce 

that ^1 ^1 where ^1 is defined in the algorithm WoMS (A2). Knowing the stopping time, we can 
easily describe the exit location since the Brownian motion is rotationnaly invariant: B^^ is then 

uniformly distributed on the sphere of radius ipAoiCi)- Hence 

(|i,||B^-J|)^i(6,x(l)), andei<T;. 

By this procedure we can construct a sequence of stopping times (^^i n > 1) and define H„ = 
^1 + . . . + H„ is the first time after such that the Brownian motion exits from a sphere 

centered in Ba ^ of radius ■i/'a^ initialized at time The moving sphere should stay in the 

domain D, so we choose 

2 ,\'^+ir(z> + 1) 



7'(^-a/Bh„_J e/(^ + l 



Using the same arguments as before and by the Markov property for the Brownian motion, we 
obtain the identities in law 

(an, n > 1) = {An, n > 1), (En, \\B^^ \\) ^ ^ = (En, xin) 
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n>l V / n>l 



Fig 1. Walk on Moving Spheres 



with En < T/ and An, x{^) defined in the algorithm WoMS (A2). Consequently defining 
= inf{n > 0;B^ ^ -C"}, the following identity yields 

(2.47) (s^^llBe^Jl) (s^.,x(iV^)) and H^, < T,. 

Step 2. Let us now estimate the difference between and T/. By (2.47) we first deduce 

(2.48) F{t) := F{ti < t) = F{Ti < t) < F%t) := P(e^. <t), t> 0. 

Furthermore, for any small a > 0, 

1 - F{t) = F{Ti >t,E^,<t-a)+ Fiji > t, H^, > t - a) 
< F{Ti > t, H^, < t - q) + P(H^, > i - a) 

(2.49) < F{Ti > t, H^, < t - a) + 1 - F'{t - a). 

At time the Brownian motion is in the e-neighborhood of the boundary dD, hence 
I — IIB- II < e. Using the strong Markov property, we obtain 

"Are 

(2.50) F{Ti > t,Ej^, <t- a) < F'{t - a) sup Fy{Ti > a). 

Since the Brownian motion is rotationnaly invariant, it suffices to choose y = {I — e,0, . . . ,0). Due 
to the convexity of dD, the following upper-bound holds 

(2.51) Fy{Ti >a)< Po( sup S^'^ < e) = Fo(2|Bi'^| < e) < 

0<t<a VZaTT 

Combining (2.48) for the upper-bound and (2.49), (2.50) and (2.51) for the lower-bound yields the 
announced estimation (2.37). □ 

2.2.2. The first time the Bessel process of index v hits a decreasing curved boundary. The algo- 
rithm developed in the previous paragraph can be adapted to the problem of hitting a deacreasing 
curved boundary. Let us define : 

(2.52) r = inf{t > : = l{t)], where I is decreasing and /(O) > 0. 
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Assumption 2.11. There exists a constant Amin > which bounds the derivative of I 



I'it) > -A^in, Vt > 0. 

The procedure then also consists in building a WoMS which reaches a neighborhood of the 
boundary. But instead of deahng with a fixed boundary as in Section 2.2.1, that is a ball of radius 
/, we shall in this section introduce the following moving boundary: the ball centered in the origin 
and of radius l{t). The arguments developed in order to prove Theorem 2.8 and Theorem 2.9 will 
be adapted to this new context. 



Algorithm (A3): 

Let us define the following positive constants 

(2.53) L = max(/(0),A„,in,^/^), ^ = y+i^^2.+2 ^(t^ + l)- 

Initialization: Set x(0) = 0, = 0, Hq = 0, Ao = Kil{0) - ^/x{0) 
The n-th step: While the following condition holds 



l{En^i) - yjx{n - 1) > e 

(denoted by C(n — 1)), we simulate Un an uniform distributed random vector on [0, 1] l-^-l"*"^, G„, a 
standard gaussian random variable and Vn a uniformly distributed random vector on . Un, Gn 
and Vn have to be independent. We then construct (.^„, H„, x(?T')) using (2.35). At the end of this 

/ \2(v+l) 

step we set An = K\l{En) - y'xin)] 

The algorithm stops when C(n — 1) is not longer satisfied: we set = and so H„ = and 
X{n) = x{n - 1). 

Outcome The exit position x(n) and the exit time. 



Let us note that the stochastic process {x{n), n > 0) is not a Markov chain since the sequence 
{An)n>o depends on both (H„,x(n)). That is why we define the following Markov chain 

Rn := {En,x{n)) eM.1 

stopped at the first time the condition C(n) is not satisfied. In the following, we shall denote 
this stopping time (number of steps of the algorithm): 



= inf{n > 0; /(E„) - < e}. 

Theorem 2.12. The number of steps of the algorithm WoMS {A3) is almost surely finite. 
Moreover, there exist a constant Cs > and eo{5) > 0, such that 

E[N^] < Csl logel, for all e < eo{6). 

Theorem 2.13. As e goes to zero, Hate converges in distribution towards t defined by (2.52) 
(with cumulative distribution function F), which is almost surely finite. Moreover, for any a > 
small enough, 

(2.54) fl - —^)F%t -a)< F{t) < F%t), for all t > 0, 

where F^(t) := F{En^ < t). 
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Proof of Theorem 2.12. 

The proof is based mainly on arguments already presented in Theorem 2.8. So we let the details of 
the proof to the reader and focus our attention to the main ideas. 

(1) The process x{f^)) is a homogeneous Markov chain and the associated operator is given by 
(2.55) Pt,.f:= I f{s,y)F({t,x),{ds,dy) 

J{s,y)ml ^ 

where / is a non-negative function and P( {t, x), (ds, dy) I is the transition probability of the chain. 



The chain starts with (Ho,x(0)) = (0,0) and is stopped the first time when /(S„) — ^/xJji) < e. 
Classical potential theory ensures that 



n=0 



is solution of the following equation 

r (Pit, x) = Pt^x<j) + g{t, x), (t, x) G D', 
\ ct){t,x) = 0, y{t,x) G dD", 



(2.56) 



with = {{t,x) G : l{t) — ^/x < e}. For the particular choice g = 1, we obtain (p{t,x) = 
Kt^x[^^]j therefore the averaged number of step is given by (/)(0, 0). 

(2) In order to point out an upper-bound for the averaged number of steps, we use a comparison 
result: we are looking for a function U{t, x) such that 



(2.57) 



J U{t,x) > Pt^xU + 1, y{t,x)eD'' 
\ U{t,x) > 0, V(t,x) G dD^. 

For such a particular function, we can deduce (j){t,x) < U{t,x). Let us define 

U{t,x) = clog ( I l{i(i)_-^>o}, 

with some constant c > which shall be specified later on. The positivity assumption on the 
boundary dD^ is trivial. Moreover since I is a decreasing function, (2.55) implies 

(2.58) Pt,.U= I U{s,y)F({t,x),ids,dy)) < [ U{t,y)F({t,x), {ds,dy) 

By using the Taylor expansion, we get 

(2.5yj U[t,y)<U[t,x) 2{l{t)-^y 3(/(t)_^)3' i^^V) ^ 

Using similar arguments and similar bounds as those presented in Lemma 2.10, the odd powers in 
the Taylor expansion don't play any role in the integral (2.58). Therefore we obtain 

Pt,.U < U{t, x)-^- I [f'^^^ Ut, x), {ds, dy) 
< U{t,x) - - / —— j=-^ ^(r)dro-i(dn) 



2 Jo J si m - V^Y 
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where is given in (2.46) and K is defined by (2.40) with A = k{1{s) — y/x)'^^'^^^\ We have now 
Pt,.U<U{t,x)-'-tpl(^-JIL_^^^ (^1^^ QV(-logr)Mr)dr). 

An appropriate choice of the constant c leads to (2.57). Finahy we get 

E[N^] < C/(0,0) = c log(/(0)/e). 

□ 

Proof of Theorem 2.13. 

The arguments are similar to those developped for Theorem 2.9, the extension of the convergence 
result to curved boundaries is straightforward. That is why we shall not repeat the proof but just 
focus our attention on the only point which is quite different. We need to prove that the Markov 
chain i?„ = (H„,x(?T')) stays in the domain = {{t,x) : < x < l^{t)} so that the hitting time 
T defined by (2.52) satisfies r > H^ye- In other words, if the Markov chain = (S„,x(?^)) for the 
n-th step is equal to {s,x), then Rn+i should belong to {{t,x) : t > s, x < /^(t)}. In the WoMS 
setting, for t > s, this means that the ball centered in x and of time-dependent radius ipAit — s) 
always belongs as time elapses to the ball centered in of radius l{t). We recall that 

A = k{1{s) - y/^f^'^+^l 

Therefore we shall prove that 

(2.60) yt > S, ijA^t -s) + ^/x< l{t). 

In fact, due to Assumption 2.11 and the definition of ipA, it suffices to obtain 

(2.61) '^l;A{t-s)<l{s)-^/^-A^in{t-s), ys<t<s + W^, 
where 

1 1 

j[ \ 2u+2 / K \ 2^+2 



W = — — = — — (l(s) -^/^) = ^ (Us) - V^). 

Vr(z. + i)2V Vr(z^ + i)2V ^ ' lV5 ' 

Due to the definition of the constant L, we have 



The right hand side of the preceding inequality is the positive root of the polynomial function 
P{X) = A^inX^ + V2(i/ + l)/eX - {l{s) - Vi). We deduce that P{W) < 0. By (2.44) and 
P{W) < 0, we obtain 



</(s)_^_A^inI^2 



<l{s)-^/x-A^in{t-s), \/s<t<S + W^. 
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Finally we have proved (2.61) and so (2.60). □ 
If Assumption 2.11 is not satisfied then it is difficult to have a general description of an iterated 
procedure in order to simulate hitting times. However the particular form of the function ipa defined 
by (2.36) permits to describe a WoMS algorithm for the square root boundaries. Let us therefore 
consider the following functions: 



(2.62) V'a(t) = ^2tlog ^^^^ "^^^^^^^ and f{t) = Vf^t, 
well defined for t < to := min ^a-^, where a = a(r(i/ + 1)2*^)"^. 

The algorithm is essentially based on the following result (the constants r and u associated with 
the hitting problem of a square root boundary for the Bessel process shall be specified in the proof 
of Proposition 2.15) 

Lemma 2.14. Let us define 

(2.63) ^'^('■'''^ = 2 (iTTTj r(i/ + l) e-"/2, r>0,u>0. 
If a = Fy{r, u) then 

(2.64) ipa{t)<fit), for all 0<t<a^. 

Proof. We are looking for a particular value a depending on both r and u such that the following 
bound holds: ipa{t) < f{t), for all < t < tQ. Since t < to, it suffices to prove that 

2tlog^ <r-ut^ g{t) := t(2log + uj < r. 

Let us compute the maximum of the function g on the interval [0, to]i with tQ fixed, 

g'{t) = 2log-^ + u-2{u + l). 

We have 

g'{t) = 0^1og^ = z. + l- -^ = aexp |- - I. - l|. 
In other words the maximum of the function g is reached for 

— r A 

tmax = exp < — — 1> 

^l2(z/ + l) / 

and is equal to 

5(imax) =5max = 2(z^ + l)a-+ie2(-+i) . 
Choosing (^max < r we obtain in particular (2.64) that means 



For ao = i r(i/ + l)e-"/2^ we get (2.64) since to = ■ □ 
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The aim is now to construct an algorithm which permits to approximate the hitting time of 
the square root boundary. Therefore we consider a Bessel process of dimension 6 which hits the 
decreasing curved boundary f{t) given by (2.62). 



Algorithm (A4) — the square root boundary: l{t) = \//3o — (3it with /3o > 0, /3i > 0. 

Let K g]0, 1[. 

Initialization: Set x(0) = 0, = 0, Sq = 0, = kF„{Po,Pi). 
The (n + l)-th step : While the following condition holds 

/(H„) — a/ x("-) > sfdenoted by C(n)), 

we define 

(2.65) An = nF,[mn) - a/^)',/3i(i - 7^)) 

where Fy is defined by (2.63), and we simulate Un+i an uniform distributed random vector on 
[0, l]'-'^-!"''^, Gn+i a standard gaussian random variable and Vn+i a uniformly distributed random 
vector onS^. Un+i, Gn+i andVn+i have to be independent. We then construct (.^„+i, x(?T'+1)) 
using (2.35). 

The algorithm stops when C(n) is not longer satisfied : we set in+\ = and so r.n+i = S„ and 
X{n + 1) = x{n). 



Proposition 2.15. The statements of Theorem 2.12 and Theorem 2.13 are true for the Algo- 
rithm (A4) associated with the square root boundary. 

Proof. All the arguments developped for decreasing boundaries with lower-bounded derivatives 
are easily adapted to the square root boundary. We let the details to the reader and focus our 
attention to the following fact: the stochastic process (H„,x(?T'), n > 0) stays in the domain 
defined by 

= {{t,x) G Rl : l{t) - > 0}. 

In the WoMS setting, for t > s, this means that for (H„,x(f^)) = is,x) G the following step 
leads to •\/x('^ + 1) < K'^n+i)- By (2.35), it suffices to prove that 

(2.66) ^/x + '^l^A{t) <l{s + t), for alH G {u > : min(/(s + u),V'A(n)) > 0}, 

with A = kF^((/(s) - V^)2,/3i(l - f^)), since x{n + I) < (TxR + V'A„(en+i))^- By Lemma 
2.14 and due to the coefficient k, we have 



^P^it)<^ilis)-V^y-^^[l-^)t. 

Hence 

{l{s + t)- V^f - ijAitf > [VKs?-Pit -V^Y- {l{s) - + /3i (1 - 

> 2^(l{s) - V/2(s)-/3ii) -^j^t> 0. 
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This leads directly to (2.66). 



□ 



Remark 2.16. The whole study points out a new efficient algorithm in order to simulate Bessel 
hitting times for given levels or curved boundaries. We can use this algorithm in two generalized 
situation: 

(1) We have assumed that the Bessel process starts from the origin. Of course the procedure 
presented here can also he applied to Bessel processes starting from x > 0. It suffices to 
change the initialization step ! 

(2) We focused our attention to the Bessel process hut we linked also the initial prohlem to the 
exit time of a 5-dimensional Brownian motion from a hall of radius I. The Algorithm (Al) 
extended to higher dimensions can also he used in order to evaluate exit times of general 
compact domains whose boundary is regular. 



3. Numerical results. In this part we wih ihustrate the pre- <" 

vious results on some numerical examples. Let us figure first an °^ 

outcome of our algorithm, the exit position from a sphere with ° 

0.2 

radius depending on time. The figure opposite is giving this result ^ 

for an radius / = 1 and a precision e = 10"'^. -o.z 

Let us compare our algorithm with existing results. Consider the -o* 

classical Euler scheme for a Brownian motion and evaluate the first ""^ 
hitting time and hitting position from a disk with given radius. 



First of all we can verify that the distribution of the hitting time for the WoMS algorithm matches 
the distribution of the hitting time of a given level for the 2— dimensional Bessel process. Figure 2 
gives this result for a starting disk with radius 1, a precision e = 10"'^ and a number of simulations 
N = 20000. In the Euler scheme the time step is At = 10~^. 




Fig 2. Distribution of the hitting time (Euler scheme and WoMS algorithm (Al)) - Histogram of the angle for the 
exit position 



We can also test the fact that the exit position is uniformly distributed on the circle. In order to 
do this we can evaluate the angle of the exit position in our WoMS procedure and show that it is 
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an uniform distributed random variable with values in [— 7r,7r]. Figure 2 also shows the histogram 
of the result for a disk of radius 1 an e = 10~^ and 20000 simulations. 

Let us now present a simulation with Algorithm (A2). We consider the hitting time of the level 
/ = 2 for the Bessel process of index = 2 and we illustrate Theorem 2.8 by Figure 3. The 
curve represents the averaged number of steps versus the precision e = 10^*^, k = 1, . . . ,7. We can 
observe that the number of steps is better than suspected since the curve is sublinear. We obtain 
the following values (for 7 = 0.9 and 100000 simulations in order to evaluate the mean). 




Fig 3. Averaged number of step of Algorithm (A2) versus e 

Finally we present the dependence of the averaged number of steps of Algorithm (A2) with 
respect to the dimension of the Bessel process. For that purpose, we simulate hitting time of the 
level I = 2 with e = 10~^, 7 = 0.9, 50 000 simulations for each estimation of the averaged value, 
and the dimension of the Bessel process takes value in the set {2, 3, . . . , 18}. 




Fig 4. Averaged number of step of Algorithm (A2) versus 5 = 2v + 2 



4. Application to the Cox-Ingersoll-Ross process. We now aim to estimate the hitting 
time of a level I > for {Xf , t > 0), a Cox-Ingersoll-Ross process. The CIR process is the solution 
of the following stochastic differential equation: 



(4.1) 



dXf = (a + bXf)dt + c\/ \Xf\dBt 
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where xq > 0, a > 0, 6 G M, c > and {Bt, t > 0) is a standard Brownian motion. We denote here 
6 = Aa/c^. 

We wiU first recall a connection between this stochastic process and {Y\t), t>0), the square of 
the Bessel process BESQ(5), solution of the equation 

(4.2) Y^{t) = yo + 5t + 2 ^* ^\Y^{s)\dBs, t > 0. 



Lemma 4.1. The CIR process has the same distribution as (Xt, t > 0) which is defined by 
(4.3) 



Xi = e''*y^(g(i-e-^*; 



, 4b ^ 



Xo = Y'{0), 

where Y is the square of a Bessel process in dimension 5 = Aa/c? (see [15]). 

Proof. Let us only sketch some ideas of the proof. Let Y^{t) be the square of the (5-dimensional 
Bessel process. By applying Ito formula, we get the stochastic differential equation satisfied by the 
process Xt 



2 

dXt = bXtdt + e''d(Y(^{l 



— (? — 

= bXt dt + bS- dt + 2e''y le-f-^X^I dB^^^_^_,,^ 

(4.4) = {a + bXt)dt + 2eTj{X^\dB^ 

where 6 = Aa/c^. Let us remark that 



^(1 - e-"*) = J%\s)ds, with Pit) = I e-f . 
We can deduce that there exists a Brownian motion {Pt, t > 0) such that 

for all t > 0. With this notation, the equation (4.4) writes 



dXt = (a + bXt) dt + 2e- J \Xt\pit) d/3, 



= ia + bXt)dt + cyJ\Xt\d(3t, 

and Xq = Y{0). This proves that the process {Xt, t > 0) has the same distribution as the CIR 
process given by (4.1). □ 

Let us consider the hitting time of a given level I for the CIR process and denote it by T/. This 
time is defined by: 

Ti = inf{s > 0;Xf = 

The previous Lemma 4.1 gives also an equivalence (in distribution) connecting the hitting time of 
the CIR process and the hitting time of the square of a (5-dimensional Bessel process. 
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Proposition 4.2. The hitting time Ti of a level I > for a CIR process has the same distri- 
bution as log ^1 — ^T^)^ where 

r^ = inf{t>0;y^(i) = /(l-^t)}, 
and is the square of a Bessel process of dimension 5 = Aa/c'^ . 

Proof. By using Lemma 4.1, has the same distribution as T; given by 

(4.5) T, = inf{s>0; y5(^(l - e"^^)) = /e"^^}. 

Define t = ^(1 — we have two situations: 

First case: If fe < let s = r]{t) where 



r/(t) = -ilog(l-^t), fori>0. 



The map rj is a strictly non-decreasing function, and we get thus 

Ti = inf {r?(t); t > 0, Y^t) = I (l - ^t^^ = t] (^inf {t > 0; Y\t) = ^ (l " ^i) }) ■ 

Second case: If 6 > 0, let also s = In this case the variable t takes its values only on the 
interval 



0,t5j.So 



,2 



The condition < t < ^ can be omitted in the estimation of the infimum as the boundary to hit : 
1 ~ 7^ is negative outside this interval and the Bessel process is always positive. Furthermore the 
function ij is also non-decreasing for b >0, and the result is thus obtained. □ 

Application of Algorithm (A4): 

An immediate consequence of Proposition 4.2 is that the hitting time Ti is related to the first time 



the Bessel process of dimension 6 = 4a/c^ reaches the curved boundary: f{t) = y Z^l — ^tj. We 

are able to apply Algorithm (A4) if Aa/c^ G N* and 6 > (the boundary is then decreasing). 
Let us denote by the number of steps of (A4) and Htv^ the approximated hitting time of the 
Bessel process associated with the particular curved boundary /. Combining Proposition 2.15 and 
Proposition 4.2 leads to 

for a small enough and t > 0. 
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APPENDIX A: SIMULATION OF RANDOM VARIABLES 
Let us introduce simulation procedures related to particular probability density functions. 

Proposition A.l. Let Z be a random variable with Gamma distribution Gamma(a, /?) i.e. 

1 



P(Z G dz) 



-z"-^e ^l{,>o}dz, a>0, /3 > 0. 



r(a)/3 

Then W = exp(— Z) has the following distribution 



W e dr) = ^,J^(-logr)--Vi/^-il[o,i](r)dr. 



In particular the stopping time defined by (2.16) has the same law as 



r(i/+i)2'^ 



"^^ e ^ . Here 



Z is a Gamma distributed random variable with parameters a = u + 2 and (3 = -jj^ . 

Proof. Let / be a non-negative function. Using suitable changes of variables, we obtain 



1 



r(a)/3° 
1 

r(a)/3° 



/(e-^)z"-^e~^dz 
/(r)(-logr)°-Vi/^-Mr. 



In order to end the proof it suffices to multiply W hj a constant and use once again a change of 
variables formula. □ 



We need to simulate Gamma distributed variables. Let us just recall some common facts. 

Proposition A. 2. (i) If a & N (so-called Erlang distributions) then the Gamma distributed 
variables Z has the same law as — (3 log{Ui ... Ua) , where {Ui)i<i<a are independent uni- 
formly distributed random variables. Hence W defined by W = exp(— Z) can be simulated by 

{UiU2...Ua)f^. 

(a) If a — 1/2 G N then Z has the same law as — /31og(C/i . . . C/[aj) + where (f/i)i<i<LaJ ^''"^ 
i.i.d. uniformly distributed random variables and N is an independent standard Gaussian r.v. 
(see, for instance, [5] chapter IX. 3). 
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